Method for generating a correlated sequence of variates with desired marginal distribution for testing a model of a communications system

ABSTRACT

To provide a time series having correlated interarrival times for use in traffic and queueing studies, an independent identically distributed random number stream is first transformed to a sequence of correlated uniform variates by a special modulo-1 autoregressive operation and then further transformed into a sequence of correlated variates with a prescribed marginal distribution from which there is developed the desired time studies.

FIELD OF INVENTION

This invention relates to a process for generating a correlated sequenceof variates, and more particularly to such generation for use inqueueing systems for modelling traffic studies of communication systems.

BACKGROUND OF THE INVENTION

For traffic studies in communication systems, it has been fairlystandard to use queueing models that assume that customer interarrivalintervals and service times are independent, each being modelled as arenewal process. These standard models are easy to simulate and undersuitable assumptions are even analytically tractable.

Unfortunately, they are very often poor models of real-life systemswhere correlations do exist, for example, bursty interval streams aretypical examples of high correlation. Additionally, video signalsinclude considerable correlation from line to line and frame to frame.In practice, there is a dramatic degradation of performance induced byincreasing correlations. These correlations when suitably accounted for,can yield improved system design. In any event, it is advantageous toinclude known correlations when modeling real systems.

Typically, simulations have used a pseudo-random number generator togenerate real sequences which model independent identically distributed(iid) uniform variates. Simulations for generating dependent sequenceswith specific correlations have received less attention.

A preferred generation method for dependent sequences of random variatesshould satisfy the following criteria:

a) Period: the generated sequence should have a reasonably long period;in particular, if the correlated sequence is generated by transformingan underlying random number sequence, then the period of the twosequences should be comparable.

b) Complexity: the time and space complexity required to implement themethod on a computer should be reasonable functions of theircounterparts in the underlying random number generator. Space complexityis rarely the issue, so one usually wants fast methods, especially whenlong sequences are required.

c) Coverage: the full range of theoretic correlations, that is (-1, 1),should be attainable by a suitable adjustment of the method'sparameters. Thus, if ρ_(min) and ρ_(max) are the extremal lag-1autocorrelations, it is required that correlation in the range (ρ_(min),ρ_(max)) be attainable. The endpoints are exempt since they may only beattained by degenerate methods.

d) Behavior: a method should generate a wide range of sample pathbehaviors as a function of its parameters.

The autocorrelation structure ρ_(n) of a real random stationary sequence{X_(n) } with finite second moment, is a measure of (linear) dependencecaptured by the sequence ##EQU1## of the lag-n correlation coefficients,where X_(n) are the random variables, and where μ_(x) is the common meanand σ_(x) ² is the common variance of the X_(n).

Various methods that have been used hitherto to generate correlatedvariates with uniform marginals can be broadly classified into threeapproaches (the goal is always to generate a stationary sequence {X_(n)} with a marginal distribution function F).

Approach 1: Generate X₀ ˜F and then continue recursively, generatingX_(n+1) from X_(n) through some autoregressive scheme which is closedunder F (i.e., preserves F). This approach is the most popular one. Itenjoys the advantage that the attendant autocorrelation function istypically very easy to derive. On the other hand, this approach is adhoc in the sense that one needs to devise a separate autoregressivescheme for each F. For some families of distributions, finding theappropriate autoregressive scheme can be quite daunting.

Approach 2: Fix a tractable marginal distribution F* (normal is atypical choice), and generate the corresponding sequence {X_(n) ^(*) },with an autoregressive scheme closed under F*. Then, the sequence{F*(X_(n) ^(*))} is uniform and the target sequence {X_(n) } is obtainedvia the inversion method, i.e., X_(n) =F⁻¹ (F*(X_(n) ^(*))). Here theautocorrelation structure has been transformed twice: once by F* andonce by F⁻¹. Nevertheless, the autocorrelation function of {X_(n) ^(*) }can be computed in principle, though numerical difficulties could bepresented by certain F. This approach is more general than the previousone, but it shifts the difficulties to the numerical calculation of thetransformed autocorrelation structure.

Approach 3: This is a special case of the previous approach when F* isalready chosen to be uniform, thereby obviating the first transformationin Approach 2. Nevertheless, the major simplication brought about bythis choice of F* merits viewing it as a separate approach.

For our novel technique, we shall follow Approach 3, and so wehenceforth focus on generation methods for correlated uniform sequences{U_(n) } and on their autocorrelation structure.

A general though unspecific procedure to generate a c.i.d. (correlatedidentically distributed) sequence of random variables {X_(n) } is totake an i.i.d. (independent identically distributed) Uniform (0,1)sequence {Z_(n) } (generated approximately on a computer via apseudo-random number generator), and transform it into a correlateduniform sequence {U_(n) }. If a correlated sequence with another(non-uniform) marginal distribution is desired, the sequence {U_(n) }may be further transformed, say by transformations D_(n) , yielding thetarget sequence {X_(n) }, where X_(n) =D_(n) (U_(n)). For historicalreasons the D_(n) shall be referred herein as distortions. Typically,D_(n) will implement an inversion method; that is D_(n) =F_(n) ⁻¹ whereF_(n) is the distribution function of the X_(n). To keep notationsimple, we shall henceforth deal with stationary sequences with a commondistortion D. Note, however, that D transforms the autocovariancestructure, i.e., it maps the autocovariance function of {U_(n) }

    Ψ.sub.U (τ)=Cov[U.sub.0, U.sub.τ ]=E[U.sub.0 U.sub.τ -1/4, τ=0, 1,                                               (1.1)

to the autocovariance function of {X_(n) }

    Ψ.sub.X (τ)=Cov[X.sub.0, X.sub.τ ]=E[D(U.sub.0)D(U.sub.τ)]-μ.sub.X.sup.2, τ=0, 1,(1.2)

where μ_(X) is the common mean of the X_(n). In particular, Ψ_(X)(0)=σ_(X) ² where σ_(X) ² is the common variance of the X_(n) (assumedfinite throughout the discussion).

In practice, constraints on a correlated sequence are stated in terms ofthe target sequence {X_(n) }. These requirements will usually includethe form of the marginal distribution of {X_(n) } as well as the firstfew autocorrelations. The need is to devise an auxiliary uniformsequence {U_(n) } and implement it on a computer such that transformingeach U_(n) to X_(n) via D would comply with the prescribed constraintson {X_(n) }. The term "controlling the autocorrelations" refers to anability to satisfy autocovariance constraints. Thus, one way to controlautocorrelations is to devise generation methods of uniform sequences{U_(n) }, derive their autocovariance function (1.1), and then proceedto find the mapping Ψ_(U) →Ψ_(X) (induced by D) and its inverse mappingΨ_(X) →Ψ_(U).

The literature contains considerable work on generating target sequences{X_(n) } directly. Such methods are tailored to a specific targetdistribution and have easily derived autocorrelation functions. However,they may require more than one transformation to conveniently generatearbitrary marginals (the first transformation being the distributionfunction and yielding a uniform variate). Since one usually needs tofollow the uniform distribution route, working with uniform sequences ismuch desired. The generation problem of correlated uniform sequences{U_(n) } was addressed by a number of workers in the prior art. Forexample, workers have proposed minification (moving minimum) andmaxification (moving maximum) processes. The present invention employs anew method for generating Markovian sequences of random variables withuniform marginals.

SUMMARY OF THE INVENTION

In particular, a feature of the present invention is a novel method forgenerating sequences of random variables with uniform marginals, to bedescribed as "Transform-Expand-Sample" (TES). This method, to bedescribed in more detail subsequently, uses as inputs an independentidentically distributed random number stream of uniform marginals and apair of specified control parameters, and subjects such inputs to amodulo 1 autoregressive scheme that yields a sequence of correlatedvariates with a uniform marginal distribution.

Often, when the sequence of variates is to be used in the modelling oftraffic studies with varying characteristics, it is advantageous thatthe sequence have a prescribed marginal distribution. In this case, byan additional computer transformation on each TES-generated variate inthe uniform marginal distributions, there can be formed a correlatedsequence with a wide range of marginal distributions.

More specifically, a random number generator is used to produce asequence of independent identically distributed variates from which oneis chosen to serve as the initial variate or seed in a modulo oneautoregressive operation. To it is added as a second component in suchoperation, an innovation variate that is constructed from anotherindependent identically-distributed variate, which may be a succeedingvariate in the sequence. In particular, the innovation variate isconstructed to have a desired probability distribution in an intervalchosen appropriately to introduce the desired autocorrelation. Thefractional part of the sum of the addition serves both as the firstvariate of the desired output sequence and as one element of asubsequent modulo one addition operation including a new innovationvariate constructed in the manner of the first innovation variate. Thisprocess is repeated recursively to derived the desired sequence ofcorrelated variates with a uniform probability distribution.

Moreover, to widen the range of autocorrelation functions available, amodified version of the just-discussed process is available in which theproduced sequence of variates is used to form a modified sequence inwhich the even-numbered variates in the sequence retain their originalvalue but the odd-numbered variates are replaced by variates whose valueis one minus their value in the original sequence.

When a desired non-uniform marginal probability distribution is desired,this derived sequence can be further modified to achieve the desiredmarginal probability distribution by use of standard distortiontechniques.

In practice, for given traffic data with a characteristicautocorrelation function, one would proceed to generate a TES sequencethat approximates the data marginal distribution and its autocorrelationfunction, and then use such a sequence in a Monte Carlo simulation. Themodelling task of finding an appropriate TES model for the given trafficdata is made easier by the availability of fast numerical methods forcalculating TES autocorrelation functions.

BRIEF DESCRIPTION OF THE DRAWING

The single FIGURE of the drawing is a flow chart of the basic operationsof the process of the invention for use in traffic studies.

DETAILED DESCRIPTION OF THE INVENTION

There is a natural affinity between Uniform (0,1) sequences and modulo-1arithmetic. In fact, uniformity is closed under modulo-1 addition in arather general way to be made precise later. This fact motivates theconsideration of modular processes, defined here as processes obtainedby applying modulo arithmetic as the final operation. It turns out thatone can define very simple modular processes leading to tractable Ψ_(u)and Ψ_(x), yet possessing interesting and varied sample path behavior,as well as favorable coverage and period properties.

The basis for our novel approach is provided by a simple fact, (Lemma 1below) establishing that uniformity is closed under modulo-1 addition.

Let x = max{integer n:n ≦x} be the integral part of x, (x real) and letx =x- x be the fractional part of x. Note that x ≧0 for all real x andthat bold brackets and braces are used to denote the integral part andfractional part operators, so as to distinguish them from theirsyntactic counterparts.

Lemma 1 (General Iterated Uniformity)

Let U˜Uniform (0,1) and X= U+V where V is independent of U but has anarbitrary distribution. Then it can be shown that X˜Uniform(0,1).

Based on this, we have devised a class of simple modular processes withuniform marginals. The algorithmic construction will be called a TES(Transform-Expand-Sample) method. Specifically, we shall be interestedin two related classes of TES methods called TES⁺ and TES⁻, giving riseto modular processes {U_(n) ⁺ }_(n=0).sup.∞ and {U _(n) ⁻ }_(n=0).sup.∞respectively. These sequences (called TES sequences) are definedrecursively over a common probability space by ##EQU2## where U₀ ⁺ =U₀ ⁻=U₀ ˜Uniform(0,1), and {V_(n) }_(n=0).sup.∞ is a stationary innovationsequence, i.e., a sequence of i.i.d variates where each V_(n) (theinnovation variate) is independent of the histories {U₀ ⁺, . . . ,U_(n-1) ⁺ } and {U₀ ⁻, . . . , U_(n-1) ⁻ } and can have a uniformprobability distribution or any other desired probability distribution.The reason for the introduction of TES⁻ is the need to provide coverageof the entire correlation coefficient range for marginally uniformsequences: TES⁺ methods (appropriately parameterized) realize all valuesof the first autocorrelation in the range [0,1], while their TES⁻counterparts cover the range [-1,0]. In a similar vein, distorted TESsequences will cover realizable positive and negative ranges of theirfirst autocorrelation.

From now on and without further mention, we shall consistently appendthe proper superscript (plus or minus) to distinguish between similarmathematical objects associated with TES⁺ and TES⁻, respectively. Forexample, ρ_(U) ⁺ (n) and ρ_(U) ⁻ (n) will denote the autocorrelationcoefficient functions of {U_(n) ⁺ } and {U_(n) ⁻ }, respectively. Themore precise equivalent notation ρ_(U+) (n) and ρ_(U-) (n) is morecumbersome and will not be used. Objects common to {U_(n) ⁺ } and {U_(n)⁻ }, or those for which the distinction is not important, will carry nosuperscript. For example, if {X_(n) ⁺ } and {X_(n) ⁻ } are thestationary distorted sequences given by X_(n) ⁺ =D(U_(n) ⁺) and X_(n) ⁻=D(U_(n) ⁻), then μ_(X) and σ_(X) ² will denote their common mean andcommon variance, in notational agreement with (1.2 ).

Using simple properties of modulo 1 addition, we can equivalentlyrepresent (1.3) as ##EQU3## is the sum of n i.i.d. variates (S₀ =0).From (1.4), the analogous representation of U_(n) ⁻ is ##EQU4## Lettingf(x) be the common density function of the innovations V_(n) anddenoting by f_(n) (x) the density function of S_(n), it is clear thatf_(n) (x) is the n-fold convolution of f. We denote by f₀ the Diracfunction δ(x) which constitutes the identify under convolution.Throughout this paper we shall consistently append a tilde to a functionsymbol h(x) to denote its Laplace transform ##EQU5## With this notation,the convolutions f_(n) (x) imply the transform relations

    fn(s)=[f(s)].sup.n, n=1,2,                                 (1.8)

The sample path relationship of Eq. (1.4) between {U_(n) ⁺ } and {U_(n)⁻ } will reduce the amount of work needed to analyze the latter, bydeducing many of its properties from the former. For example, Eq. (1.4)implies that

{U_(n) ⁻ }_(n=0).sup.∞ ={U₀, 1-U₁ ⁺, U₂ ⁺, 1-U₃ ⁺, . . . } leadingimmediately to the important simple relations

    ρ.sub.U.sup.- (τ)=(-1).sup.τ ρ.sub.U.sup.+ (τ), τ=0, 1,                                               (1.9)

    Ψ.sub.U.sup.- (τ)=(-1).sup.τ Ψ.sub.U.sup.+ (τ), τ=0, 1,                                               (1.10)

As mentioned before, although Lemma 1 allows the innovations {V_(n) } tohave any distribution, it is clear that computational and analyticalconsiderations will steer us to "simple" innovations. I shall focusmostly on

    V.sub.n ˜Uniform(-L,R), n=1, 2,                      (1.11)

where the pairs (L,R) satisfy 0≦L, R≦1, L+R≦1. Since each pair (L,R)completely specifies a TES method, we shall denote by TES⁺ (L,R) andTES⁻ (L,R) the corresponding algorithmic construction of the resultant{U_(n) ⁺ } and {U_(n) ⁻ } respectively. In the sequel, we shall oftenfind it convenient to switch from the (L,R) parameterization of TESmethods to the equivalent parameterization (α, φ) given by ##EQU6## Inthis case we write TES⁺ (α, φ) and TES⁻ (α, φ).

The simplicity of {V_(n) }_(n=0).sup.∞ in (1.11) is evident from thefact that it can be easily constructed as a linear transformation of ani.i.d. Uniform(0,1) sequence {Z_(n) }_(n=1).sup.∞. We identify {Z_(n) }with the values of a pseudo random number stream. Specifically,

    V.sub.n =-L+αZ.sub.n, n=1,2,                         (1.14)

so the Laplace transforms f_(n) (S) for the sums S_(n) in (1.6) become##EQU7## or equivalently, under the (α,φ) parameterization ##EQU8## Weremark again that the special choices in (1.14) and (1.15) will beassumed only for practical calculations and only where explicitlystated.

TES methods have a simple geometric interpretation as random walks onthe unit circle. For example, TES⁺ (L,R) constructs U_(n+1) ⁺ from U_(n)⁺ by forming a "circular interval" (arc) around U_(n) ⁺ of length L toits "left" and length R to its "right" and then samples uniformly forU_(n+1) ⁺ in this arc. The TES⁻ (L,R) construction is similar exceptthat the arc is constructed about 1-U_(n) ⁻. In the alternativeparameterization (α, φ) of (1.12)-(1.13), α is the total length of thearc whereas φ is a measure of the offset of the arc from symmetricstraddle. Symmetric straddle corresponds to φ=0 (i.e., R=L), so in TES⁺(α, 0) the variate U_(n+1) ⁺ is generated with equal probability 1/2 tothe left or to the right of U_(n) ⁺, and the sequence {U_(n) ⁺ } haszero net drift around the unit circle. For φ>0 (i.e., R>L), U_(n+1) ⁺ ismore likely to lie to the right of U_(n) ⁺ than to its left, and thesequence {U_(n) ⁺ } has a net clockwise drift around the unit circle.The case φ<0, will analogously produce a net counter-clockwise drift.

TES methods enjoy the following properties.

a) A TES method presupposes the availability of a standard i.i.d.uniform random number stream {W_(n) }. However, its simplicity makes itsuitable for computer generation, and the complexity of generating thecorrelated sequence {U_(n) }, equals that of generating {W_(n) } plus asmall constant. Furthermore, its geometrical interpretation (to bediscussed later) shows that the period of {U_(n) } is at least as longas that of {W_(n) }.

b) For any prescribed -1≦r≦1, there is a TES method with Corr[U_(n),U_(n+1) ]=r, n≧1. Furthermore, TES methods give rise to autocorrelationfunctions ρ_(n) whose magnitude in n is both monotonic (for φ=0) andoscillatory (for φ≠0); in this sense, TES methods improve on theminification/maxification methods whose resultant autocorrelationfunctions are strictly monotonic.

c) While a basic TES method generates a stationary process with uniformmarginals, it can be easily modified to generate non-stationary uniformsequences {U_(n) } with arbitrary step dependent correlations r_(n)=Corr[U_(n), U_(n+1) ]. In fact, the r_(n) could be modulated by time orby some other stochastic mechanism.

d) TES methods are particularly suitable for simulating random cyclicphenomena. Random resets occur when the point 0 is crossed in eitherdirection.

A more detailed discussion of TES follows.

This discussion will be organized as follows. First, there is presentednotational conventions and other preliminaries. TES methods are formallyintroduced next. Then, we compute and investigate analytically the firstautocorrelation of TES generated uniform sequences, whereas higher lagautocorrelations are investigated empirically by simulation. Finally, wetouch on model fitting issues with TES sequences. The autocorrelationfunction of TES sequences and their transformations (via the inversionmethod) can be derived through Fourier methods; numerical calculationsof autocorrelation functions of TES sequences were found to be inexcellent agreement with their simulation counterparts including thoseto be discussed. This shows that TES induced autocorrelation functionscan be computed numerically from convergent series representations.

In the following discussion, the following notation is used. Theindicator function of a set A is denoted 1_(A) (·); however, when A isan event, we usually omit the (sample point) argument. If x is a realnumber, then x denotes its integral part and x its fractional part(equivalently its value modulo 1); notice the bold type used todistinguish the bracket and brace operators from their ordinarysyntactic meaning. Specifically,

     x =sup {integer n: n≦x}

     x =x- x

so x is always non-negative (even for negative x).

It is convenient to view the interval [0,1) as a circle, obtained byjoining the points 0 and 1 together, say, at the bottom of the circle.The topology of the circular unit interval provides a convenientgeometric model for the wraparound effect of modulo 1 arithmetic. To seethat, observe that for any two points x, y ε [0,1) there are twodistance-like functions

    d.sub.l (x,y)= x-y  (left-distance)                        (2.1)

    d.sub.r (x,y)= y-x  (right-distance)                       (2.2)

We note in passing the relation d_(l) (x,y)=d_(r) (y,x), which for x≠yis also a direct consequence of the identity

     x + -x =1, xε(0,1).                               (2.3)

Neither d_(l) (·) nor d_(r) (·) is a metric since they are notsymmetric. However, a metric d_(c) can be defined on [0,1) in a naturalway by

    d.sub.c (x,y)=min(d.sub.l (x,y), d.sub.r (x,y)).           (2.4)

We shall refer to d_(c) as the circular distance on [0,1). It is alsouseful to define relations ≦_(c) and <_(c) on [0,1) (called the circularpositioning relation and strict circular positioning relation) by

    x≦.sub.c y iff d.sub.c (x,y)=d.sub.l (x,y) iff d.sub.c (x,y)=d.sub.r (y,x)                                                     (2.5)

    x<.sub.c y iff x≦.sub.c y and x≠y             (2.6)

where the left-hand side of (2.5) reads "x is to the left of y"(equivalently, "y is to the right of x"), and similarly for (2.6). Whileevery two points in [0,1) are comparable under the circular positioningrelation, it does not constitute a total ordering due to the absence oftransitivity. Nevertheless, it does provide a sense of localdirectionality compatible with the usual geometric interpretation ofclockwise rotation as movement to the right and counter-clockwiserotation as movement to the left. Furthermore, crossing point 0 ineither direction represents the modulo 1 operation, or more pictorially,wraparound action. A wraparound interval a,b on the circular unitinterval with endpoints a and b in [0,1), a≠b, is ##EQU9## where boldbrackets and parentheses are used to distinguish wraparound intervalsfrom ordinary ones. In words, a,b consists of all points on the circularunit interval from a and clockwise to b (including a but excluding b),possibly wrapping across point O. It is convenient to define 0,1 =[0,1)and somewhat inconsistently a,a ={a} (singleton set). The definition ofother circular intervals with open/closed endpoints is similar.

We begin with a general description of TES (Transform-Expand-Sample)methods. Subsequently, we progressively specialize to TES methods ofparticular interest. A general TES method operates within the followingmathematical setup. Let,

    θ:[0,1)→[0,1) be a bijection (1--1 and onto mapping),

    δ:[0,1)→2.sup.[0,1) be a set-valued function from [0,1) to subsets of [0,1).

Denote,

    D.sub.y =δ(y)

    C.sub.x ={yε[0,1):xεD.sub.y }.

Thus, C_(x) is the set of all y ε[0,1) such that D_(y) "covers" theprescribed x ε[0,1). Letting λ(·) denote the Lebesgue measure on [0,1]we make a number of mild regularity assumptions (all measurabilitystatements are with respect to the Lebesgue measure):

3.a) θ(·) is measurable.

3.b) D_(y) is measurable for almost every y.

3.c) C_(x) is measurable for almost every x.

Definition 1: A general TES method iterates a three step loop togenerate U_(n=1) given θ(·), δ(·) and U_(n) =y:

Step 1: Transform y→θ(y).

Step 2: Expand θ(y)→D.sub.θ(y) =δ(θ(y)).

Step 3: Sample uniformly in D.sub.θ(y) to obtain U_(n+1), so that theconditional density of U_(n+1) given U₀, . . . , U_(n) depends only onU_(n) via the set D.sub.θ(U.sbsb.n.sub.), viz. ##EQU10##

Our interest in TES methods is motivated by the previously mentionedLemma of Iterated Uniformity, now discussed in more detail.

Let U˜Uniform(0,1) and let V be obtained by applying a TES method as inDefinition 1. Suppose further, that for some αε(0,1) we have:

a) λ(D_(y))≡α for almost every yε[0,1).

b) λ(C_(x))≡α for almost every xε[0,1).

Then V˜Uniform(0,1).

Proof: Directly from the construction, V has a conditional density(given U), ##EQU11##

Therefore, the (unconditional) density of V is ##EQU12## as required.

The reason for the term "Iterated Uniformity", is motivated by thefollowing.

Corollary 1: The sequence {U_(n) }_(n=0).sup.∞ (obtained by theiterative application of a TES method on an initial U₀ ˜Uniform(0,1)) isstationary with Uniform[0,1) marginals. Furthermore, it forms a MarkovProcess (discrete parameter, continuous state) with transition densitygiven by (3.1), which is also the joint density of each pair (U_(n),U_(n+1)).

Lemma 1 can be generalized, by dropping conditions a) and b) provided westill have ##EQU13## The problem is then to find suitable θ(·) and δ(·)such that the attendant D.sub.θ(y) still satisfies (3.2).

We now proceed to specialize TES to an important subclass characterizedby sets D_(y) of the wraparound interval type. We, henceforth, focus onTES methods determined by triples TES(θ, L, R) where

a) θ(·) is a measurable bijection as before.

b) L, R ε[0,1] are constants such that L+R=αε[0,1].

The sets D.sub.θ(y) are constructed as D.sub.θ(y) = θ(y)-L , θ(y)+R ,i.e., as wraparound intervals of, respectively, length L to the left andlength R to the right of θ(y), yielding a wraparound interval D.sub.θ(y)of total length α. Note that as long as θ(y)-L≧0 and θ(y)+R≦1, the setD.sub.θ(y) forms one contiguous interval in the linear view of [0,1); ifeither θ(y)-L<0 or θ(y)+R>1, wraparound will yield a D.sub.θ(y) which isthe union of two possibly non-contiguous intervals anchored at thepoints 1 and 0 respectively. In the circular view of [0,1), the setD.sub.θ(y) always forms a contiguous interval.

It is convenient to interchange the L and R parameters with anequivalent parameter φ defined as ##EQU14## Recalling that,

    α=R+L,                                               (3.4)

the equivalence of (α,φ) and (L,R) follows immediately from (3.3) and(3.4) justifying the interchangeable notation TES(θ,L,R) or TES(θ,α,φ).Observe that φ is a measure of how θ(y) is offset from the center of theinterval D.sub.θ(y) (alternatively, φ measures how D.sub.θ(y) is rotatedrelative to the inverval straddling θ(y) symmetrically). Strictlyspeaking, the use of φ as a structural parameter characterizing a TESmethod is redundant, since φ can be absorbed in θ(·). However, it can bejustified as a matter of convenience since it simplifies the explanationof the qualitative dynamics of the random walk {U_(n) }. To see that,observe that φ can be represented as ##EQU15## Thus, p and q are theprobabilities of sampling U_(n+1) to the left and to the right ofθ(U_(n)), respectively. Their difference φ then represents the drift (or"velocity") of {U_(n) } as it meanders around the circular unitinterval. Positive φ correspond to positive velocity (right or clockwisedrift) while negative φ correspond to negative velocity (left orcounter-clockwise drift). The interpretation for φ=0 is obvious. Theextremal cases φ=±1 correspond to maximal velocities (φ as defined in(3.3) is standardized by division by α).

In the remainder of the discussion we focus on two fundamental mappingsθ denoted θ⁺ and θ⁻ and defined by

    θ.sup.+ (y)=y  (identity)                            (3.6)

    θ.sup.- (y)=1-y  (reflection)                        (3.7)

Definition 2: TES⁺ (L,R) (or TES⁺ (α,φ)) is the TES method TES(φ⁺, L,R)giving rise to a sequence U⁺ ={U_(n) ⁺ } defined recursively as follows:

a) U₀ ⁺ =U₀ ˜Uniform(0,1).

b) Given U_(n) ⁺, the next iterate U_(n+1) ⁺ is sampled uniformly on

    D.sub.U.sbsb.n.spsb.+ =  U.sub.n.sup.+ -L , U.sub.n.sup.+ +R  ,

i.e., ##EQU16##

Definition 3: TES⁻ (L,R) (or TES⁻ (α,φ)) consists of alternating TESmethods TES(θ⁻, L, R) and TES(θ⁻, R, L) giving rise to U⁻ ={U_(n) ⁻ }defined recursively as follows:

a) U₀ ⁻ =U₀ ˜Uniform(0,1).

b) Given U_(n) ⁻, the next iterate U_(n+1) ⁻ is sampled uniformly on##EQU17##

The attendant correlation coefficients of lag n are similarly denoted(n>0)

    Corr[U.sub.0.sup.+, U.sub.n.sup.+ ]=ρ.sub.n.sup.+,

    Corr[U.sub.0.sup.-, U.sub.n.sup.- ]=ρ.sub.n.sup.-.

For n=1 we simply write ρ⁺ and ρ⁻, omitting the subscript.

The reason for emphasizing the TES⁺ and TES⁻ methods is that togetherthey generate every lag -1 correlation; in fact, we show later that thecorrelation coefficients ρ⁺ (α,φ) and ρ⁻ (α,φ) as a function of α and φmap (α,φ) ε[0,1]×[-1,1] onto the intervals [0,1] and [-1,0],respectively. First we establish that these are indeed TES methodsgiving rise to uniform marginals.

Proposition 1: TES⁺ and TES⁻ are TES methods satisfying Lemma 1.

Proof: Both θ⁺ and θ⁻ are obviously Lebesgue measurable. Since^(D).sub.θ.spsb.+.sub.(y) and ^(D).sub.θ.spsb.-.sub.(y) are unions of atmost two intervals, they are clearly measurable. Directly fromdefinition, they have Lebesque measure α. Finally for TES⁺, C_(x) ⁺ =x-R , x+L ; and for TES⁻, C_(x) ⁻ = -x-A , -x+B , where A=R and B=L onodd iterates, whereas A=L and B=R on even iterates. Since both C_(x) ⁺and C_(x) ⁻ are wraparound intervals, they are measurable. Each has

Lebesgue measure α because they can be represented as translations ofthe intervals D.sub.θ.spsb.+.sub.(x) and D.sub.θ.spsb.-.sub.(x), on thecircular unit interval.

The proof of the last part of Proposition 1 becomes trivial whenconsidering the geometrical interpretation of wraparound as described inSection 2. It is then easy to see that the wraparound intervals C_(x) ⁺and C_(x) ⁻ are (rigid) translations of the basic interval [0,α).

Now we calculate the lag-1 autocorrelation of TES⁺ and TES⁻. We begin byshowing that it is enough to investigate TES⁺ since the relation between{U_(n) ⁺ } (induced by TES⁺) and {U_(n) ³¹ } (induced by TES⁻) isantithetic in the following sense.

Proposition 2: There is a construction of {U_(n) ⁺ } and {U_(n) ⁻ } overa common probability space such that ##EQU18##

Proof: It is enough to prove for n=0,1,2 as the rest follows by simpleinduction.

For n=0, U₀ ⁻ =U₀ ⁺ =U₀ by definition.

For n=1, let W₁ ˜Uniform (0,α) be independent of U₀. Directly, fromdefinition,

    U.sub.1.sup.+ = U.sub.0 -L+W.sub.1  .

Define W₁ '=α-W₁ ; then W₁ '˜Uniform(0,α) independent of U₀.Furthermore, with probability 1, ##EQU19## Finally, for n=2, let W₂˜Uniform(0,α) be independent of U₀, U₁ and W₁, Then,

    U.sub.2.sup.- =U.sub.1.sup.+ -L+W.sub.2 =U.sub.2.sup.+,

as required.

Corollary 2: ρ_(n) ⁻ =(-1)^(n) ρ_(n) ⁺, n=1,2, . . . , since byProposition 2 ##EQU20##

Corollary 3: It is straightforward to conclude from Proposition 2 therepresentations ##EQU21## where U₀ ˜Uniform(0,1) and the W_(i)˜Uniform(0,α) are i.i.d. and independent of U_(j), 0≦j≦i.

We now proceed to derive the lag-1 autocorrelation ρ⁺.

Theorem 1: Let {U_(n) ⁺ } be generated by TES⁺ (α,φ). Then, ##EQU22##

Proof: Let α be restricted initially to the interval (0,1). We use therepresentation

    U.sub.0.sup.+ =U.sub.0,U.sub.1.sup.+ =U.sub.0 -L+W.sub.1

where U₀ ˜Uniform(0,1) and W₁ ˜Uniform(0,α) are independent. Now,##EQU23##

Next, evaluating the partial expectations we get ##EQU24##

Routine algebraic manipulations now yield the result for αε (0,1). Thelimiting cases ρ⁺ (0,φ)≡1 and ρ⁺ (1,φ)≡0 are correctly captured by(4.1), since for α=0, D_(x) ={x}, and for α=1, D_(x) =[0,1).

Corollary 4: For any αε[0,1], φε[-1,1], we have ρ⁺ (α,φ)=ρ⁺ (α,-φ),since ρ⁺ (α,φ) is a function of φ². Thus, TES(θ⁺, L, R) and TES(θ⁺, R,L) yield the same lag-1 autocorrelation.

Theorem 2: Let {U_(n) ⁻ } be generated by TES⁻ (α,φ). Then, ##EQU25##

Proof: Immediate from the relations

    ρ.sup.- (α,φ)=-ρ.sup.+ (α,-φ)=-ρ.sup.+ (α,φ),

where the first equality is justified by Corollary 2, and the second oneby Corollary 4.

The analytical behavior of ρ⁺ (α,φ) and ρ⁻ (α,φ) is straightforward.Since they differ only by opposing signs, we concentrate on ρ⁺ (α,φ).Solving the equation ρ⁺ (α,φ)=0 for α with φ viewed as a parameter, onefinds that its zeros α₁ and α₂, are real non-negative for all φ. Infact, ##EQU26## Notice that α₁ ≧1 implies φ² ≦1/3. For φ² >1/3, ρ⁺ (α,φ)is negative for ##EQU27## reaching a minimum at ##EQU28## Now, the pairof curves ##EQU29## partition the open (α,φ) rectangular R=(0,1)×(-1,1)into sets F, G and H defined by ##EQU30## such that ρ⁺ (α,φ) is positiveon H, negative on G and vanishes on F. For ρ⁻ (α,φ), the roles of G andH are reversed.

From the viewpoint of simulation applications, this section can besummarized as follows.

Corollary 5: Let {r_(n) } be any sequence of reals such that -1<r_(n)<1, and let U₀ ˜Uniform(0,1) be given. Then, for each n≧1, there is aTES method T_(n) that operates on a uniform variate U_(n-1) yielding auniform variate U_(n) such that

    Corr[U.sub.n-1,U.sub.n ]=r.sub.n, n≧1.

Furthermore, if 0<r_(n) <1, then T_(n) =TES⁺ (α,φ) for any pair (α,φ) εH satisfying the equation ρ⁺ (α,φ)=r_(n), i.e., ##EQU31## or T_(n) =TES⁻(α,φ) for any pair (α,φ) ε G satisfying the equation ρ⁻ (α,φ)=r_(n),i.e., ##EQU32## Similarly, if -1<r_(n) <0, then T_(n) =TES⁻ (α,φ) forany pair (α,φ) ε H satisfying ρ⁻ (α,φ)=r_(n), or T_(n) =TES⁺ (α,φ) forany pair (α,φ) ε G satisfying ρ⁺ (α,φ)=r_(n). In either case (-1<r_(n)<1), for each -1≦φ≦1 there is a unique solution for α such that (α,φ)lies in the region H.

The correlation coefficients r_(n) need not be constant, because a TESmethod generates U_(n+1) as a function of U_(n), r_(n) and the valueW_(n) from an i.i.d. Uniform(0,1) random number stream. Consequently,the sequence {r_(n) } need not even be deterministic; one can simulatemodels where a random mechanism generates the correlations r_(n),thereby adding another modulation component to the {U_(n) }.

The computation of higher lag autocorrelations is more complicated. ForTES⁺ sequences, these representations yield for n=1, 2, . . . ##EQU33##where B_(n) (X) is the Bernoulli function of order n obtained from theBernoulli polynomials by their periodic extension from the interval[0,1) to the entire real line (see Knopp Theory and Application ofInfinite Series, published by Blackil and Son, London and Glasgow(1946)). The ρ⁻ (α,φ) can be easily obtained from the last equation viaCorollary 2.

Finally, some aspects concerning the use of TES methods as models forcorrelated data sequences will be discussed, also touching briefly onsome of their statistical properties as gleaned from computer-generatedruns.

First, it is to be noted that restricting the sequence {U_(n) } to haveuniform marginals is usually without loss of control over the attendantlag-1 autocorrelation. More specifically, let {X_(n) } be a realstationary random sequence with marginal distribution F. Suppose thateach X_(n) =F⁻¹ (U_(n)) is obtained by an inversion method from anunderlying sequence {U_(n) } with Uniform(0,1) marginals (see Bratley,et al A Guide to Simulation, published by Springer-Verlag (1987),Devroye Non-Uniform Random Variate Generation, published bySpringer-Verlag (1986)). For a prescribed -1≦ρX≦1, require that

    Corr[X.sub.n, X.sub.n+1 ]=ρX, n≧0,

and find (α,φ) such that ρ(α,φ) for {U_(n) } gets transformed to ρX for{X_(n) }. Now, it is known that the extremal correlations are attainedas follows:

    Corr[F.sup.-1 (U.sub.n), F.sup.-1 (U.sub.n)] (positive extremal)

    Corr[F.sup.-1 (U.sub.n), F.sup.-1 (1-U.sub.n)] (negative extremal)

Since we can find a TES method to generate {U_(n) } with any lag-1autocorrelation -1≦ρ≦1, we expect that any lag-1 autocorrelation canalso be induced on the {X_(n) } spanning the interval between theirextremal values (ρ_(min), ρ_(max)). Note that this interval may beproperly contained in (-1,1); for example, the extremal negativecorrelation for the exponential case is about 0.645. Observe that ρ⁺(α,φ) and ρ⁻ (α, φ) are monotonic in (α,φ) in the region H (but not inG). Consequently, ρX (α,φ) is also monotonic in (α,φ) in the region Happroaching its extremal values where ρ⁺ (α,φ) or ρ⁻ (α,φ) approachtheirs. The transformation ρX (α,φ)→ρ(α,φ) is then monotonic in (α,φ), afact useful in generating conversion tables for each F.

Recall that Corollary 5 implies that a TES method has two "degrees offreedom", namely, α and φ. Since Eqs. (4.1) and (4.2) for ρ⁺ (α,φ) andρ⁻ (α,φ) are quadratic in both α and φ, their behavior is readilyunderstood. Clearly, the magnitude of ρ⁺ and ρ⁻ decreases in themagnitude of α and φ in the region H. However, the α parameter provideshere the primary controlling effect; in other words, for fixed φ, anylag-1 autocorrelation can be obtained by an appropriate choice of α. Theeffect of φ on the correlation is secondary. If we fix α and set theextremal values φ=0 and |φ|=1 in ρ⁺ (α,φ) or ρ⁻ (α,φ) we get ##EQU34##and this difference attains its maximal value at α=1/2 resulting in amaximal difference of 3/8. However, φ controls an important aspect ofthe sample path behavior of {U_(n) }, namely, its drift (or "velocity")around the circular unit interval. Recall that in Eq. (3.5), p and q canbe taken to measure the left (counter-clockwise) and right (clockwise)drift of {U_(n) } so that φ=q-p can be interpreted as the net drift ofU_(n+1) relative to U_(n) around the circular unit interval.Intuitively, φ measures how often {U_(n) } crosses a given point in agiven direction. Note that the function θ(·) itself can provide a driftcomponent; however, for θ⁺ (·) and θ⁻ (·) (i.e., identify andreflection), this component is zero and φ can be justifiably argued tofully account for the (standardized) velocity of {U_(n) }.

The sample paths of TES⁺ are cyclical in nature. Qualitatively speaking,they exhibit a discontinuity in the neighborhood of the point 0 due towraparound. When crossing 0 clockwise on the circular interval, onepasses from relatively large values to relatively small ones, andconversely when crossing counter-clockwise. This effect in the samplepaths will make sense when modeling physical phenomena with a cyclicalnature where "resetting" punctuates runs of uptrending or downtrendingobservations, though the runs need not be monotonic. Naturally, this isnot always justifiable from a modeling viewpoint. A possible extensionof these principles is to find TES methods whose sample paths look more"continuous" than those introduced here. This will mean that theunderlying process {U_(n) } is prevented from crossing point 0 bytreating it as a reflecting boundary. Mathematically, this will requireother solutions for θ(·) and δ(·) in Eq. (3.2). It is tempting to try alinear autoregressive scheme (without the modulo 1 operation). Morespecifically, consider

    U.sub.n+1 =U.sub.n -a(U.sub.n)+b(U.sub.n)W.sub.n+1         (5.1)

where 0≦a(y),b(y)<1,yε[0,1], are continuous deterministic functions andW_(n+1) ˜Uniform(0,1) is the innovation, subject to the constraints

    y-a(y)≧0, yε[0,1]                           (5.2)

    y-a(y)+b(y)<1, yε[0,1]                             (5.3)

(constraints (5.2)-(5.3) ensure that D.sub.θ(y) never straddles thepoint 0). However, the scheme (5.1) with the constraints (5.2)-(5.3) canonly be satisfied trivially, as the following Proposition makes precise.

Proposition 3: If the sequence {U_(n) } in (5.1) is to have Uniform(0,1)marginals, then necessarily a(y)=y and b(y)=1 for almost every yε[0,1).

Proof: For 0≦y≦1,0≦x<1, ##EQU35## By assumption,

    P{U.sub.n+1 ≦x}=x, xε[0,1).

Equating the righthand sides above, we get ##EQU36## For this functionalequation to hold we must have ##EQU37## From (5.2) and (5.4) we concludethat a(y)-y=0 or equivalently a(y)=y, almost everywhere. From (5.2) and(5.3) it follows that b(y)≦1 which together with (5.5) implies b(y)-1=0or equivalently b(y)=1, almost everywhere.

Thus, the autogressive scheme (5.1) reduces to U_(n+1) =W_(n+1) whichcorresponds to the trivial case of the i.i.d. Uniform(0,1) stream.

We now illustrate how a fitting methodology might be employed to fit aTES method to observed data. The discussion is kept at an intuitivelevel and is meant to outline the issues rather than to treat them indepth. Suppose we have a sample sequence {U_(n) }_(n=0) ^(N) assumed tobe stationary and with values in [0,1). We would like to capture thelag-1 autocorrelation as well as its drift through a TES method T. Let rbe an estimate of the lag-1 autocorrelation, and let φ=q-p be anestimate of the drift φ. Then, use T=TES⁺ or T=TES⁻ according as r isnon-negative or non-positive, respectively, and estimate α by solving(for α) as explained in Corollary 5. This could be satisfactory if thedata satisfy

    d.sub.c (U.sub.n, U.sub.n+1)≦α, n≦0.   (5.6)

If not, one could consider a modified TES method defined as aprobabilistic mixture ##EQU38## where 0<π<1 and V_(n+1) ˜Uniform(0,1) isindependent of U₀, . . . ,U_(n). Note that Corr[U_(n),U_(n+1) ]=π·ρ(α,φ)where ρ(α,φ) is the correlation coefficient in the constituent pure TESmethod. The effect of the mixture is to introduce regeneration pointswith probability 1-π and, therefore, to dilute the magnitude of the pureTES lag-1 autocorrelation by a factor π; the resulting {U_(n) } consistof independent runs of pure TES variates of geometrical length (withmean ##EQU39## In general, a mixed TES method gets around the degeneracyaspect of pure TES methods which stipulate (5.6) with probability 1.

Estimating φ can sometimes be difficult. A simple approach is toestimate p and q by the relative frequency of local increases ordecreases, i.e., ##EQU40## where # is the cardinality of the followingset. This can be expected to work well when the data exhibit high lag-1autocorrelation so that successive values are relatively close in thesense of the metric d_(c). Otherwise, it is hard to distinguish betweenincreases (decreases) and the effect of the modulo 1 operation resultingfrom crossing the point 0 in either direction. Thus, for relativelylarge α (α≧1/2) it is not clear how to assign a pair of successiveobservations into p or q, even when the simple estimators p and q aboveare replaced by the more realistic ones ##EQU41## where ≦_(c) is thecircular positioning relation defined in Eqs. (2.5)-(2.6).

If there is to be provided a prescribed autocorrelation in the TESgenerated sequence, it is important to be able to measure the desiredcorrelation to be prescribed.

We now proceed to describe fast numerical methods for calculating theautocorrelation (and autocovariance) functions of TES sequences, as wellas their transformed variety.

A basic result we shall employ to study autocorrelation properties isthe Poisson Summation Formula which states that if g is a differentiablefunction such that both ##EQU42## exist, then ##EQU43## (i=√-1). Therighthand side of (7.1) is the Fourier series of the periodic function##EQU44## and the integrals are the corresponding Fourier coefficients.Note that (7.1) can be written in the form ##EQU45## Finally, thefollowing simple proposition will be frequently used.

Proposition 1

Let h(x) be a real function whose Laplace Transform ##EQU46## exists.

Then, h(ir) and h(-ir) are complex conjugate for any real r.

Proof

The complex exponential pairs e^(ir) =cos(r)+i sin(r) and e^(-ir)=cos(r)-i sin(r) are obviously complex conjugates. Since h(ir) andh(-ir) are obtained by integrating these with respect to h(x)dx, itfollows that h(ir) and h(-ir) are also complex conjugates.

Now we derive the transition structure of general TES sequences. Denoteby g.sub.τ⁺ (y|x) the conditional density of U.sub.τ⁺ given U₀ =x, andby g.sub.τ⁻ (y|x), the conditional density of U.sub.τ⁻ given U₀ =x. Inthe next two lemmas we proceed to calculate the conditional densitiesabove; note that we make no assumptions on the f_(n) (x).

Lemma 2 {U_(n) ⁺ } is a stationary Markov process with transitiondensity g.sub.τ⁺ (y|x), where for τ=0, 1, . . . ##EQU47## Proof {U_(n) ⁺} is stationary Markovian by construction; see Eq. (3.1).

For τ=0, both sides of Eq. (3.1) reduce to the Dirac function δ(y-x).For τ>0 we can write by virtue of (1.5), ##EQU48## Since the functione^(-s{x+)σ} is periodic in σ (period 1), we can rewrite the above as##EQU49## where the interchange of summation and integration in thesecond equality is justified by uniform convergence, and in the thirdequality because a Fourier Series can be integrated termwise; note alsothat we used the Poisson Summation Formula (7.2) in the last equality.Since 0≦x<1, we can decompose the last integral as follows ##EQU50##whence, after adding the integrals, ##EQU51##

But ##EQU52## is the Laplace Transform of the function e(y)=e^(i2)τvy1(y≧0), so inverting c.sub.τ (s|x) yields ##EQU53## The effect of e^(-s)in d.sub.τ (s|x) is translation by 1, so inverting d.sub.τ (s|x) yields##EQU54## by the periodicity of the complex exponential. Eq. (3.1) nowfollows by calculating g.sub.τ⁺ (y|x)=c.sub.τ (y|x)-d.sub.τ (y|x) fromthe expressions above. Finally, Eq. (8.2) follows from Eq. (8.1) bynoting that the term for v=0 evaluates to 1; furthermore, termscorresponding to index pairs (-v,v) are complex conjugates byProposition 1, so summing them yields twice their real part.

Lemma 3

{U_(n) ⁻ } is a stationary Markov process with transition densityg.sub.τ⁻ (y|x), where for τ=0,1, . . . ##EQU55## or in real form,##EQU56##

Proof

{U_(n) ⁻ } is stationary Markovian by construction; see Eq. (8.4).Directly from Eq. (8.4) we have ##EQU57## Since the case for τ even wasalready computed in Lemma 2, it remains to prove (3.2) for the case τodd. Substituting 1-y for y in (8.1) now yields Eq. (8.3) immediately##EQU58## where the last equality follows by the periodicity of thecomplex exponential. Eq. (8.4) follows from Eq. (8.3) by an argumentanalogous to that at the end of Lemma 2.

Lemmas 2 and 3 easily specialize to TES processes with uniformmarginals.

Corollary 1

For TES⁺ (α,φ) and TES⁻ (α,φ), the corresponding transition densitiesare for τ≧0, ##EQU59##

Proof

From Eq. (1.16), we have ##EQU60## where the last equality follows fromEuler's Formula for the sine function. The rest follows by substituting(8.7) into (8.2) and (8.4), respectively.

It is of interest to compare TES processes to their time-reversedcounterparts which are also stationary Markovian. Let {U_(n) ⁺}_(n=0).sup.τ be a TES⁺ subsequence and define its time-reversed version{U_(n) ⁺ }_(n=0).sup.τ by U_(n) ⁺ =U.sub.τ-n⁺. Now from Eq. (1.5)

    U.sub.τ.sup.+ ={U.sub.0 +S.sub.τ }.                (8.8)

Using modulo 1 arithmetic rules we can isolate U₀ as U₀ ={U.sub.τ⁺-S.sub.τ }, or equivalently in terms of the time-reversed subsequence

    U.sub.τ.sup.+ ={U.sub.0.sup.+ -S.sub.n }.              (8.9)

A comparison of (3.8) and (3.9) shows that a time-reversed TES⁺ sequenceis a TES⁺ sequence with an innovation sequence {V_(n) } related to{V_(n) } by V_(n) =-V.sub.τ-n. We refer to the V_(n) as reverseinnovations. For the special case of uniform innovations, reverseinnovations become

    V.sub.n =L-αZ.sub.n =α-R-αZ.sub.n =-R+α(1-Z.sub.n).(8.10)

Since {1-Z_(n) } is an i.i.d. Uniform(0,1) sequence, the effect ofreverse innovations is to change TES⁺ (L,R) to TES⁺ (R,L); equivalently,TES⁺ (α,φ) under time reversal becomes TES⁺ (α,-φ) by virtue of (1.13).

Consider now the time-reversed subsequence {U_(n) ⁻ }_(n=0).sup.τdefined by U_(n) ⁻ =U.sub.τ-n⁻. For even τ, the situation is analoguousto TES⁺ since the even index TES⁻ subsequence is an even index TES⁺subsequence. For odd τ, the situation is quite different. In fact, fromEq. (1.7)

    U.sub.τ.sup.- = -U.sub.0 -S.sub.τ                  (8.11)

so that U₀ = -U.sub.τ⁻ S.sub.τ , or equivalently in terms of thetime-reversed subsequence

    U.sub.τ.sup.- = -U.sub.0 -S.sub.τ                  (8.12)

Clearly, the representations (8.11) and (8.12) are identical. Therefore,we reach the somewhat surprising conclusion that the odd indexsubsequence of a TES⁻ process is time-reversible in the sense that thetransition density from U₀ to U.sub.τ⁻ is identical to the transitiondensity from U₉₆ ⁻ to U₀. Notice that for TES⁻ V_(n) =V.sub.τ-n, SO{V_(n) } has the same probability law as {V_(n) } (permuting the indicesof an i.i.d. sequence has no effect on the joint distribution).

To check our conclusions, we compare g.sub.τ⁺ (y|x) or Eq. (8.1) and g⁻(y,x) of Eq. (8.3) to their time-reversed counterparts g.sub.τ⁺ (y|x)and g.sub.τ⁻ (y|x). For TES⁺, the effect of reverse innovations is toreverse the sign of the argument in density transforms; that is, weshould substitute f.sub.τ (-i2πv) for f.sub.τ (i2πv) in Eq. (8.1). Wethen have the relation ##EQU61## To see that, note that reversing y andx in the sum in (8.13) yields the same sum as in (8.1) (except that theorder of summation is changed). For TES⁻ and τ odd, we have immediately

    g.sub.τ.sup.- (y|x)=g.sub.τ.sup.- (y|x)=g.sub.τ.sup.- (x|y)           (8.14)

as expected.

Consider the sequences {X_(n) ⁺ }_(n=0) .sup.∞ and {X_(n) ⁻ }_(n=0).sup.∞ where X_(n) ⁺ =D(U_(n) ⁺) and X_(n) ⁻ =D(U_(n) ⁻) are bothassumed to have finite second moment. Now we derive the autocovariancefunctions of {X_(n) ⁺ } and {X_(n) ⁻ }. We first obtain theautocovariance function Ψ_(x) ⁺ (τ) of the sequence {X_(n) ⁺ } with theaid of Lemma 2.

It can be shown that for τ=0,1, . . . ##EQU62##

It can also be shown that for τ=0,1, . . . ##EQU63##

An interesting consequence of Theorems 1 and 2 is that replacing theinnovation sequence {V_(n) } by its reverse counterpart {V_(n) }, whereV_(n) =-V.sub.τ-n, leaves the autocorrelation structure unchanged. Tosee that note that the effect of reverse innovations is to replacef.sub.τ (i2πv) by f.sub.τ (-i2πv) everywhere in Theorems 1 and 2. ButEq. (9.2) and (9.4) do not change because Re[f.sub.τ (-i2πv)]=Re[f.sub.τ(i2πv)] (from Proposition 1 we take the real part of complexconjugates). In particular, for uniform innovations, Eq. (8.9) impliesthat reverse innovations change TES⁺ (L,R) to TES⁺ (R,L). But Eq. (1.13)implies that this is equivalent to replacing φby -φ. We conclude thatTES⁺ (α,φ) and TES⁺ (α,-φ) have the same autocovariance function Ψ_(x) ⁺(τ ); similarly, TES⁻ (α,φ) and TES⁻ (α,-φ) have the sameautocorrelation function Ψ_(x) ⁻ (τ). The transition structure (and thesample paths) are, of course, quite different, except for odd index TES⁻subsequences.

Finally, we point out that Theorem 1 and 2 can be modified toaccommodate non-stationary TES sequences. If the innovations {V_(n) }are non-stationary, then Eqs. (9.1)-(9.4) will still hold. If, however,we replace D by a a distortion sequence {D_(n) }, then Eqs. (9.1) and(9.2) must be modified. Specifically, one must replace |D(i2πv)|² by D₀(i2πv)D.sub.τ (-i2πv) and D(i2πv)² by D₀ (i2πv)D.sub.τ (i2πv).Furthermore, the real form Eqs. (4.2) and (4.4) can no longer bededuced, since their proof depends heavily on stationarity.

The general formulas in Theorem 1 and Theorem 2 can be specialized intwo ways: by specializing the density f of the innovations {V_(n) }, orby specializing the distortion D to yield arbitrary marginals. Weproceed to specialize Theorem 1 and 2 through a sequence of corollaries.

Uniform Innovation and General Distortion

We begin with the important special case of uniform innovations.Computationally, this is the most useful specialization of Eqs.(9.1)-(9.4). Recall that the uniform innovation assumption implies thatthe V_(n) are of the from (1.14) and the corresponding f.sub.τ (s) aregiven by (1.16) and (8.7). Define

    D(i2πv)=a.sub.v +ib.sub.v,v≧1                    (10.1)

Then for each v≧1,

    |D(i2πv)|.sup.2 =a.sub.v.sup.2 +b.sub.v.sup.2(10.2)

    and

    Re[D(i2πv).sup.2 ]=a.sub.v.sup.2 -b.sub.v.sup.2.        (10.3)

Corollary 2 ##EQU64## Corollary 3 ##EQU65## Proof

It is to be noted that Eqs. (10.4) and (10.5) provide the basis forpractical numerical calculations of Ψ_(x) ⁺ (τ) and Ψ_(X) ⁻ (τ), for TESmethods with uniform innovations.

The case of uniform innovation and linear distortion D(x)=cx+d (c, dreal) corresponds to Uniform(c,d) marginals. A simple calculation yields##EQU66##

Corollary 4

For the uniform innovation and linear distortion case ##EQU67## Thespecial case of the identity distortion D(x)=x, obtained for c=1 andd=0, yields the Uniform(0,1) TES sequences {U_(n) ⁺ } and {U_(n) ⁻ }.Now, we use Eqs. (9.2) and (9.4) to give explicit expressions for theentire autocorrelation coefficient function ρ_(U) ⁺ (τ); theautocorrelation coefficient function ρ_(U) ⁻ (τ) is easily deduced withthe aid of Eq. (1.9).

The uniform innovation assumption in (1.14) means that each f_(n) (x)has Laplace transform given by (1.15) or (1.16). The representations ofΨ_(U) ⁺ (τ) and Ψ_(U) ⁻ (τ) are summarized in the following corollarieswith the aid of two auxiliary functions ##EQU68##

Corollary 5

The representations for Ψ_(U) ⁺ (τ) and Ψ_(U) ⁻ (τ) depend on the valueof τ (mod 4) as follows: ##EQU69## The result now follows from(10.9)-(10.10) by noting in succession that ##EQU70##

Formulas (10.11) and (10.12) in Corollary 5 reveal an interestingconnection with Bernoulli polynomials. Extending the Bernoullipolynomials, (denoted B_(n) (x)) periodically (with period 1) from theirrestriction to the interval [0,1) to the real line yields the Bernoullifunctions (denoted B_(n) (x)). It is known that the latter have therepresentation, ##EQU71##

Corollary 6

The representations of Ψ_(U) ⁺ (τ) in (10.16) and Ψ_(U) ⁻ (τ) in (10.17)can be simplified as follows: ##EQU72##

Corollary 7 ##EQU73## Proof

From σ_(U) ² =1/12.

Now we keep the setup of Section 5.1, but we specialize to thelogarithmic distortion ##EQU74## (where the ln function denotes thenatural logarithm) corresponding to the exponential distributionfunction F(x)=1-e⁻λx. If {U_(n) } is marginally Uniform(0,1), then{X_(n) }={D(U_(n))} is marginally Exponential(λ). Our calculations arebased on the identities ##EQU75## where γ≈0.57721566 is Euler's constantand ##EQU76##

are respectively the sine integral and the cosine integral. Numericaltables for these functions are available, e.g., in Jahnke and Ende.Tables of Functions with Formulae and Curves, Dover Publications (1945).Our next step is to derive |D(i2πv)|² and D(i2πv)².

Now, for v≧1, ##EQU77## Decomposing the above into real and imaginaryparts gives ##EQU78## Substituting (10.19)-(10.20) into the above yields##EQU79##

Corollary 8

For the uniform innovation and logarithmic distortion case, ##EQU80##Next we take D to be a step function with a countable range. Forexample, this situation arises naturally when an empirical histogram(frequency table) is used to estimate the distribution of a continuousor discrete variate. Formally, let N be a countable set of consecutivepositive integers starting from 1. Let {p_(k) : kεN} be probabilities,(i.e., p_(k) ≧0, kεN, and ##EQU81## Define C₀ =0 and ##EQU82## Thus,{C_(k) }_(k)εN is the cumulative distribution of {p_(k) }_(k)εN. Letfurther {r_(k) : kεN} be a sequence of reals with the interpretationthat each X_(n) ⁺ or X_(n) ⁻ assumes the value r_(k) with probabilityp_(k). Consequently, the corresponding distortion is ##EQU83## where theintervals I(k)=[C_(k-1), C_(k)), kεN partition the interval [0,1). Now,from Eq. (10.24) ##EQU84##

Corollary 9

For the uniform innovation and step function distortion case, ##EQU85##

We now illustrate the use of Corollary 9 by applying it to the case ofgeometric marginals; that is, we make the following identifications:N={integer k:k ≧1}, r_(k) =k, p_(k) =(1-q)q^(k-1) and C_(k) =1-q^(k) forsome 0<q<1, for all kεN. Substituting the above into (10.23) and (10.24)yields ##EQU86## These can be further simplified with the aid of Abel'ssummation by parts formula, which states that for a real sequence {A_(k)} such that ##EQU87## we have ##EQU88##

Corollary 10

For the uniform innovation and geometric marginal case, we have (v≧1)##EQU89##

Next we take D to be the linear distortion D(x)=cx+d but leave theinnovations arbitrary. Recall that by Lemma 1 the generated sequences{X_(n) ⁺ } and {X_(n) ⁻ } are still marginally uniform.

Corollary 11

For τ=0,1, . . . ##EQU90##

Proof

We conclude by observing that Eq. (10.30) and (10.31) are in agreementwith Eq. (1.10).

Summarizing, Theorems 1 and 2 provide the theoretical basis fornumerical computation of autocorrelation functions of transformed TESsequences. To check the efficacy of our methods, we compared themagainst Monte Carlo simulations for some of the special cases developedin above and found excellent agreement.

Finally, the basic operations of the process that has been describedabove are condensed in the flow chart that forms the single figure ofthe drawing.

As seen in the figure, there is first generated an i.i.d. random numberstream by any of the techniques known for such purpose. This stream isthen transformed to a sequence of correlated uniform variates by TESmethods involving a modulo 1 autoregressive scheme, as described. Thenthis sequence typically is transformed in known fashion to a sequence ofvariates with a desired marginal distribution. The resulting sequence isthen used as a queueing model in traffic studies in the usual fashion.

While the invention has been described with specific reference toproviding a sequence of variates with desired marginal distribution andcorrelations for traffic studies, it should be evident that the basicprinciples should be applicable to providing such a sequence forcomparable use in other applications.

What is claimed is:
 1. A computer-implemented method for testing a modelof a communications system by applying a simulated traffic patternhaving a desired correlation and detecting a response to suchapplication, the method comprises the steps:producing a sequence ofindependent identically-distributed variable signal values of uniformmarginals and choosing therefrom first and second independent signalvalues, constructing from said second chosen signal value an innovationsignal value having a prescribed marginal distribution in an interval ofsignal values chosen appropriately to achieve a desired correlation inthe simulated traffic pattern, combining the first chosen signal valueand the innovation signal value in a modulo-one addition and recoveringthe fractional value of the addition as a variate in a sequence ofautocorrelated signal values, constructing a new innovation signal valuehaving the prescribed marginal distribution from a new independentvariable signal value chosen from the sequence ofidentically-distributed variable signal values, combining the newinnovation signal value in a new modulo-one addition with the recoveredfractional value of the previous modulo-one addition and recovering thefractional value of such addition as a subsequent variate in thesequence of autocorrelated signal values, repeating in recursive fashionthe choosing, constructing, combining and recovering steps to formsucceeding variates of the sequence of autocorrelated signal values,applying the sequence of autocorrelated signal values for simulating thetraffic patterns to the model of the communications system to be tested,and detecting the output response of the model to said applied sequencefor ascertaining the real-time performance that can be anticipated ofthe communication system under the traffic pattern being simulated. 2.The computer-implemented method of claim 1 in which a random numbergenerator produces the initial sequence of independentidentically-distributed variable signal values.
 3. Thecomputer-implemented method of claim 1 further comprising the step,occurring before the applying step, of modifying by a distortionoperation the variates of the sequence of autocorrelated signal valuesrecovered from the modulo-one additions to obtain a desired marginaldistribution.